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(57) Abstract 

Method and device for filtering a digital signal with the aid of an Mm order finite impulse response filtering with a fractional delay 
d which is a fraction of the sampling period associated with the digital signal, in order to obtain an interpolated value. This involves 
predetermining (M+l)/2 window weights w* of a window function which is applied to the infinite impulse response of the ideal interpolator, 
the window function being a symmetric stepped window function. Then simplified filter coefficients (hj)wo are determined as a function 
of the fractional delay d, and the finite impulse response filtering is applied to the digital signal. Finally, the weighted result of the finite 
impulse response filtering is multiplied by a common normalisation factor (F(d))wO- Each filter coefficient determination requires just three 
operations. 
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1 

Method and device for filtering a digital signal with fractional delay 

The present invention relates to a method and device for filtering a digital signal 
in order to reconstruct a value of a signal sample at any time between two samples of 
5 the digital signal. More specifically, the invention relates to the filtering of a digital 
signal with the aid of an Mth order (M preferably odd) finite impulse response filtering 
(FIR) with a fractional delay d which is a fraction of the sampling period associated 
with the digital signal, in order to obtain an interpolated value, comprising the step of 
predetermining (M+l)/2 window weights Wk of a window function which is applied to 

10 the infinite impulse response of the ideal interpolator. 

Such filters are used in various fields where it is necessary for signals which are 
limited in bandwidth and are sampled at twice the Nyquist frequency to be recon- 
structed at any time in respect to the sampling times. These interpolation filters are also 
referred to as fractional delay filters or variable phase delay filters. 

1 5 Fields of application include, but are not limited to, synchronisation in digital 

modems, phase-controlled array antennae (for example "null-steering" antennae), asyn- 
chronous sample frequency conversion of previously sampled signals ("upsampling" 
and "downsampling"), speech coding, signal pattern recognition and temporal localisa- 
tion of patterns (for example speech recognition), periodical or non-periodical signal 

20 pattern generation in any phase desired, physical modelling of musical instruments (for 
example vibrating strings), but also the tracking of signals in digital navigation receiv- 
ers (GPS, LORAN-C, EUROFDQ, the time-continuous tuning of numerically con- 
trolled oscillators (NCO) and adaptive averaging devices (in LORAN-C), the direct 
determination of signal characteristics such as zero crossings and maxima/minima in, 

25 for example, seismic measurement data, the change in the aspect ratio (3/4 <-> 9/16) in 
the case of video processing, simulating multipath phenomena in the case of antenna 
signals in mobile receivers, the rapid determination of (local) maxima and minima in 
sampled data (for example in GPS data processing), time-discrete Hilbert transformers 
for software radio receivers, signal pattern recognition for temporal and/or positional 

30 localisation in, for example, the medical sector (ECG (heart rhythm analysis), EEG, 
etc) and the synchronisation of many asynchronous digital data streams into a single 
data stream (for example the synchronisation of compass data, gyroscope data, level 
data, roll data, pitch data, (GPS) position data, speed data, acceleration data, and data 
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regarding the change in acceleration into one navigational solution for aeroplanes), or 
the synchronisation of various asynchronous process data streams into one data stream 
for digital control systems having a plurality of inputs in all forms of industry (i.e. a 
simplification of Kalman filtering). 

5 Such filters are used as an alternative for known mathematical interpolation 

methods for determining the best curve through a particular set of equidistant points. If 
the original signal, after filtering by an ideal low pass filter has been sampled at twice 
the Nyquist frequency, the original signal can be reconstructed perfectly by convoluting 
the sample values with the impulse responses associated with the ideal low pass filter. 

10 The impulse response of the ideal low pass filter is a time-continuous sine function 
where the period of the side lobes is equal to the inverse of twice the Nyquist fre- 
quency. For the ideal case, even one convolution operation will already take an in- 
finitely long time for such a filter of infinite order M. In practice, however, the only 
possibility' is to work with causal impulse responses and a limited order of the inter- 

1 5 polation filter, the so-called finite impulse response filters ( FIR). 

The frequency response of a fractional delay filter, implemented as a finite 
impulse response filter, is generally given by 

H(z) = h 0 + h, * z 4 + h 2 * z" 2 + + h M * z" M = f> 3 * z-> 

j-o 

20 where hj represents the impulse response weights and M is the filter order. 

Finite impulse response filters with fractional delay are disclosed, for example, 
by the article "Splitting the Unit Delay" by T.I. Laakso et al., in IEEE Signal Process- 
ing Magazine, Volume 13, No. 1, pp. 30-60, January 1996. 

This publication discloses, for example, a design for a finite impulse response 
25 filter with fractional delay, where the infinite sine impulse response is multiplied by a 
bell-shaped window function. The zero values outside the bell shape ensure that the 
ultimate impulse response will be constrained or finite. The computation of the filter 
coefficients for such a finite impulse response filter is complex (for each coefficient 
determination, a value of a computation-intensive bell-shaped window function must be 
30 multiplied by a value of a computation-intensive sine function), and the error intro- 
duced by the filter is not minimised, moreover. 

Another way, disclosed in this publication, of designing finite impulse response 
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filters with3ractional delay is the so-called Lagrange interpolation. Here, the filter 
ficients hj are given by 




j = {0,1 V ..,M} 



5 

Here, too, many computation-intensive operations are required for each coeffi- 
cient, thus increasing the total operation time per interpolation of the filter. 

Since the known filter designs require many, often complex computations for 
each filter coefficient of the finite impulse response filter, it is often not possible to 
10 achieve the accuracy, speed of adjustment to new values of the fractional delay d, 
bandwidth and/or maximum allowable total filter delay which are required for a par- 
ticular application. 

It is an object of the present invention to provide a method for obtaining an in- 
terpolated value of a digital signal with a fractional delay d as a variable component of 
15 a delay of M/2 sampling times, with the aid of Mth order finite impulse response fil- 
tering, which requires the smallest possible number of operations per filter coefficient, 
so that highly accurate, rapidly adjustable filtering with fractional delay can take place 
at the highest possible bandwidth even for small values of the total filter delay M/2+d. 
This object is achieved in a method of the type defined in the preamble, charac- 
20 terised in that the window function is a symmetric stepped window function and in that 
the method further comprises the following steps: 

- determining simplified filter coefficients (hj)wo as a function of the fractional delay d 
according to the formula 



fh(M) 



25 



v 



wO 



' h (k,d) 



= (-D 



(k-i)+d 



k = {^,...,2,l} = -^-i 
|d|<i 



wO 



and subsequently applying the finite impulse response filtering to the digital signal; 
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- multiplying the weighted result of the finite impulse response filtering by a common 
normalisation factor (F(d))wo, given by the formula: 



(F(d)) w0 = 



M+i r 





f \ 




r \ 






h(M) 


+ 


h (k,d) 






O- 1 J 


w0 







which is a function of the simplified filter coefficients, in order to obtain the inter- 
polated signal value from the digital signal with fractional delay d. 

This method has the advantage that each filter coefficient globally requires only 
three operations, namely an addition, a subtraction and a division. This allows the 
10 fabrication of finite impulse response filters which can be very rapidly adjusted to the 
fractional delay d and consequently have a small total filter operation time of the size 

of 4M operation times. 

In one embodiment, the window weights of the symmetric stepped window 
function are chosen such that the resulting finite impulse response filter exhibits a fre- 
1 5 quency response having a ripple-shaped deviation of the normalised level (equiripple 
filter). In a further embodiment, the window weights of the symmetric stepped window 
function are chosen such that the resulting frequency response is maximally flat 
(Lagrange filter). This means that a choice can be made from a number of classes of 
finite impulse response filters by choosing the window weights without this having 
20 consequences in terms of the way in which the filter coefficients are calculated. 

As a result of, in one embodiment of the method according to the invention, an 
additional common factor (%-d 2 ) being added when the simplified filter coefficients 
(hj)^ and the common normalisation factor (F(d)U are determined, the method can 
also be applied to a fractional delay d equal to ±'/ 2 . In the embodiment mentioned ear- 
25 Her this value was excluded to prevent a division by zero from determining the simpli- 
fied filter coefficients. To this end, in fact, two critical divisions are replaced by two 
multiplications, the total number of operations increases slightly. 

In one embodiment of the method according to the present invention, two finite 
impulse response filtering operations are used to effect filtering which allows all fre- 
30 quencies to pass ("all pass filter"), the simplified filter coefficients being determined 
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solely fronf the window weights and the desired fractional delay d. The transfer func- 
tion of the all pass filter is given by 

h 0 +h, *z~ l + + h M *z~ M 

H -*-« W = h 0 *z" M h l *z- (M - ,) + + h M 

Hence it can be deduced that both the multiplication by the common normalisa- 
tion factor (F(d))wo and the addition of the factor (%-d 2 ) can be dispensed with for this 
filter. 

The above-described embodiments of the method according to the invention 
have the drawback that one of the operations in determining the simplified filter coeffi- 
cients and the common normalisation factor is a division. A division is an operation 
which, both when carried out by software and carried out by hardware, generally re- 
quires much more operation time than a multiplication operation. 

In a further embodiment of the method according to the invention, the division 
in the course of determining the filter coefficients and the common normalisation factor 
is replaced by an approximation with the aid of one or more multiplications. Thus, the 
finite impulse response filtering with fractional delay d can be carried out more rapidly, 
compared with the previously described embodiment, with an acceptable 
approximation error. 

In a further embodiment of the method according to the invention, the 
approximation for a main lobe of the impulse response and the approximation for side 
lobes of the impulse response take place separately in such a way that the resulting im- 
pulse response retains the exact values for the extreme values of the fractional delay 
(i.e. for d=0 and d=± l A). 

In a first approximation, the finite impulse response filtering with a fractional 
delay d, where the simplified filter coefficients and the common normalisation factor 
are calculated with the aid of multiplications, can be described by a third-order poly- 
nomial (P0=3) of the fractional delay, the maximum approximation error in the impulse 
response being 41 0" 3 . The complexity of this embodiment can, based on the approxi- 
mation by a third order polynomial, be reduced in a further approximation to an em- 
bodiment using a second-order polynomial (P0=2) with a maximum approximation 
error remaining almost the same at 6- 10°. In an alternative approximation, the filtering 
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can be described by a fifth-order polynomial (P0=5), the maximum approximation error 
being 2-l(T s . This, however, is at the cost of doubling the complexity of the filtering. 

In a further embodiment of the method according to the invention use is made 
of the inverse of the finite impulse response filtering with fractional delay d as 
described above, to determine those values of the fractional delay at which the original 
signal has zero crossings and/or extremes. Not only the point in time and the corre- 
sponding amplitude, but also the conditions in which these phenomena occur can be 
determined accurately if the order of the polynomial is less than or equal to 3 (P0<3). In 
the case of a third-order polynomial approximation, zero crossings can be determined 
via a solution of a third-order equation and extreme values according to a simpler 
solution of a second-order equation If the inverse of the second-order polynomial 
approximation of the fractional delay filter according to the present invention is used, 
zero crossings can be determined by solving a second-order equation, and extremes can 
be determined by solving a linear equation. 

A second aspect of the present invention relates to a device for implementing 
the method according to any one of Claims 1 to 7 inclusive. To this end, the device is 
provided with delay elements, addition elements, multiplication elements and inverter 
elements for filtering a digital signal with a fractional delay d which is a fraction of the 
sampling period associated with the digital signal, 

predetermined window weights w K of a window function which is applied to the infinite 
impulse response of an ideal interpolator being implemented in multiplication elements, 
simplified filter coefficients (hj)wo being implemented as a combination of addition 
elements, multiplication elements and inverter elements and 

a common normalisation factor (F(&)U being implemented as a multiplication ele- 
ment, 

in order to obtain the interpolated signal value from the digital signal with fractional 
delay d. 

As this device is capable of hardware implementation, much more rapid opera- 
tion of the filter can be obtained, compared with a software implementation. 

In a further embodiment of the device according to the present invention, the 
device is further provided with computing means in which the inverse of the finite im- 
pulse response filtering with fractional delay d is employed to determine the values of 
the fractional delay d at which the original signal has zero crossings and/or extremes. 
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A possible application of the method for filtering a digital signal in order to re- 
construct a value of a sampled signal at any time between two samples of the digital 
signal is a numerically controlled oscillator (NCO). Known NCOs (for example the 
type HSP45116 from Harris Corporation) comprise a K-bit phase accumulator and an 
electronic table containing 2 K numbers having a word size of 16 bits. In the case of a 
20-bit phase accumulator this leads to a table of 16 Mbit. With the aid of the phase 
accumulator, any frequency can be chosen in the range from 0 to f s /2 in steps of 2" K *fi. 

A further aspect of the present invention relates to a numerically controlled 
oscillator for generating digitised band-limited periodic signals, for example sinusoidal 
signals, comprising a K-bit phase accumulator and an addressable table characterised in 
that the addressable table comprises twice the oversampling (=2*b) of groups of P0+1 
polynomial coefficients which are addressed by the 2 log(2*b) most significant bits of 
the K-bit phase accumulator, 

the addressed P0+1 polynomial coefficients being used for a POth order approximation 
of a filter with a fractional delay d, formed by P0 multiplication elements and P0 addi- 
tion elements, the fractional delay d being determined by the 2 log(2*b) least significant 
bits of the K-bit phase accumulator. Since only 2*b groups of P0+1 polynomial coeffi- 
cients need be stored in the table, a much smaller table is sufficient. In the case of, for 
example, a sine signal oversampled by a factor of 16 (2b=32) and a third-order poly- 
nomial approximation (P0=3) of a 5 th order (M=5) fractional delay filter, only 128 
numbers need be stored in the table. Simulations have shown that for a second-order 
(P0=2) polynomial approximation about 1000, and for a first-order (P0=1) polynomial 
approximation about 8000 numbers should be stored for comparable quality. 

The present invention will now be explained in more detail with reference to 
examples and the appended drawings, in which 

Figure 1 shows a finite impulse response of a seventh-order filter, together with 
a window function which is applied to the infinite impulse response of the ideal low 
pass filter; 

Figure 2 shows a general schematic for the fractional delay filter, implemented 
using a finite impulse response filter; 

Figure 3 shows an example of a schematic of a fifth-order fractional delay filter, 
in which the range of the fractional delay d has been extended by the value |d| = '/i; 

Figure 4 shows an example of a schematic of a third-order fractional delay filter 
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which passes all frequencies to a same degree; 

Figure 5 gives the general schematic for a ninth order filter with 10 filter coeffi- 
cients according to an embodiment of the present invention; 

Figure 6 shows the schematic for a design of an integrated circuit for the first 
5 approximation (third-order polynomial) of the filter shown in Figure 5 ; 

Figure 7 shows the schematic for a design of an integrated circuit for the second 
approximation (fifth-order polynomial) of the filter shown in Figure 5; 

Figure 8 shows the schematic for a design of an integrated circuit of the filter 
shown in Figure 5 according to a second-order polynomial approximation, derived 
1 0 from the embodiment using a third-order polynomial approximation; 

Figure 9 shows the schematic for a design of an integrated circuit for determin- 
ing times of zero crossings and extremes and the corresponding amplitudes with the aid 
of the second-order polynomial approximation showing Figure 8; 

Figure 10 shows die schematic for a numerically controlled oscillator according 
15 to a further aspect of the present invention. 

The frequency response of a fractional delay filter, which has been implemented 
as a finite impulse response, can generally be described as 

H(z) = h 0 + h, * z- + h 2 * z- 2 + + h M * z- M = fX * z~ j 

j=0 



20 



Figure 1 shows a finite impulse response of a seventh-order filter (a fragment of 
an infinite sine function), together with the limiting window function which, in accor- 
dance with the method according to the invention, is applied to the infinite impulse 
response associated with the ideal interpolation filter. The symmetric window weights 
25 are given by w,...w 4 , the symmetry being round the centre of gravity of the impulse 
response. The filter coefficients h 0 ...h 7 are indicated by the vertical lines at times 0, T, 
7T, where T is the sampling period. The total filter delay lime is indicated in 
Figure 1 by the phase delay D and the fractional delay by the letter d. 

To filter a digital signal with the aid of the method according to the invention, 
30 the symmetric stepped window function is first applied to the infinite impulse response 
of the ideal low pass filter to produce the Mth-order interpolation filter. Because of the 
symmetry of the window function it is possible to use the parameters i and k, thus 



WO 99/60701 



PCT/NL99/00302 



allowing the filter coefficients to be divided into two equal parts j={i; M-i}. For the 
example in Figure 1, the values of t(s), j, k and i are shown below the graph. Conse- 
quently, the filter coefficients hj are derived by calculating the simplified filter coef- 
ficients (hj) w o according to the equation 



'h(k,d^ 

H 

V 

r h (k,d^ 



= (-1)' 



k-l*. 



(k-i)+d 



wO 



= (-1) 



(k-j)-d 



wO 



i = {o,i,...,M 

k = {*f,...,2,l} = -^-i 
|d|<i 



whereupon the total result of the digital signals filtered therewith are multiplied by the 
common normalisation factor (FCd))^, which is calculated using the equation 



10 



(F(d)L = 



M+l 







f \ 




fh(k,d) 


+ 


h (k,d) 






wO 







The window weights w k are determined as a function of the desired design type of the 
filter. For an equiripple filter, the window weights are given by 



15 



(2k -1) 



20 



where the coefficients (h k )mr are deterrnined via existing design software for equiripple 
filters. 

For a maximally flat filter or Lagrange filter, the window weights arc given by 
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Shown below as an example are the window weights which were calculated 
with the aid of a computer program for a finite impulse response filter where the phase 
delay ripple remains below about 0.1% and where the frequency response ripple 
remains below a value of 0.7%, based on the normalised frequency response. M is the 
filter order, and N=M+1 is the number of filter coefficients. 
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The method according to the invention is also suitable for implementing a 
0 Lagrange filter. A Lagrange filter has a maximally flat filter response, but a smaller 
frequency range compared with the above-mentioned equiripple filter. The window 
weights in this case are given by the following table: 





M 


Wl 


W 2 


W3 


w 4 


w 5 
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1/14 
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1 5 Figure 2 shows a general schematic for the fractional delay filter, implemented 

by means of a finite impulse response filter. The schematic incorporates delay ele- 
ments, indicated by Z\ multiplication elements, indicated by t>. and a summation 
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element, indicated by +. The filter coefficients (hj)^ and the common normalisation 
factor (F(d))wo are calculated in the above-described manner. Once the filter order and 
the window weights have been chosen, the filter coefficients depend solely on the frac- 
tional delay d. 

5 Figure 3 shows a schematic of an example of a fifth-order functional delay filter 

where the range of the fractional delay d has been extended by the value |d| = Y*- This 
is possible by introducing an additional common factor (Vi-d 2 ) into the computation of 
the filter coefficients. The total number of operations is increased slightly, because the 
implementation in practice of twice this additional factor in total requires four more 
10 multiplications plus two more additions. It should be noted in this context, that at the 
same time two critical divisions (critical because of division by zero) have been 
replaced by two safe multiplications. 

Figure 4 shows an example of a third-order fractional delay filter which allows 
all frequencies to pass. These filters are used for specific applications requiring a com- 
1 5 pletely flat frequency response. The filter coefficients are calculated in the same way as 
for the previously described filters, whilst the multiplication by the common normali- 
sation factor (F(d))wo can be dispensed with, as can the above-mentioned provision to 
avoid a division by zero. This last-mentioned advantage is a result of the fact that, 
compared with FIR filters, the fractional delay in all pass filters is increased by a factor 
20 of two, thus allowing the range of d to be halved and the critical division to be reliably 
avoided. Consequently, the third-order filter shown as an example in Figure 4 requires 
only three delay elements, eight multiplication elements, five addition elements and one 
subtraction element. 

According to a further embodiment of the method according to the invention, 
25 the divisions carried out to calculate the simplified filter coefficients and the common 
normalisation factor are replaced by approximations using multiplications. In a further 
embodiment of the method according to the invention, the approximations for the side 
lobes of the impulse response (k > 1) and the approximations for the main lobes (k = 1) 
of the impulse response are performed separately. Hereinafter, the filter coefficients for 
30 the approximation will be indicated by h. 

A first approximation (indicated by subscript 3 because of the third-order poly- 
nomial approximation) of the filter coefficients of the side lobes of the impulse 
response is given by 
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fh (k,d)l = (-If ^(k-i)? d)*(|-d 2 )*< 



fork>2 



where subscript i represents the rank of the coefficients from the "left-hand" term of the 
5 impulse response to which only the "upper" part of the + , i.e. a minus sign, applies and 
where subscript k represents the rank of the coefficients from the "right-hand" term of 
the impulse response to which only the "bottom" part of the T , i.e. a plus sign, applies 
and where also 



10 



C k = 



(k-i) J *F M (0) 



is a constant and the common normalisation factor F M (0) is given by 



F M (0) = w,+i*« 



M+l 
kk=2 



= W J__^ + ^_^! + ... + (-l) 
13 5 7 



M 



15 
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A second approximation (indicated by subscript 5 because of the fifth-order 
polynomial approximation) of the side lobes of the impulse response can be given by 
adding to the first approximation of the filter coefficients a term (1-B u *d 2 ), the calcula- 
tion of the filter coefficients consequently being given by 

fh (k,d)l =(-l) k - , *((k-i) : Fd)*(i-d J )*(l-Bu*J 2 )*C ll 

\j;M-i Ji 



25 



Approximations can also be derived for the main lobes, resulting in a first 
approximation of the filter coefficients given by 

(h k „ t (d,M)^ =(i+d)*(l + (i-d 2 )*CC,) 
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where only the "upper" part of the + , i.e. a minus sign, applies to the coefficients from 
the left-hand term of the impulse response and where only the "bottom" part of the + , 
i:e. a plus sign, applies to the coefficient from the right-hand term of the impulse 
5 response and where also CCi = Ci - 4. 

A second approximation is provided by the addition of the factor (1-Bi*d ), 

resulting in the formula 

(h w (d,M)l -(* TdH + & -d 2 )*(l-B, *d 2 )*CC,) 

10 

The correction factors B k can be calculated almost exactly for a specific value of the 
fractional delay, namely d = 0.361. At this value of the fractional delay, the first 
approximation error reaches a maximum as the filter order M increases. 

Again, two situations can be distinguished, viz. the correction factors B k (k > 1) 
1 5 for the side lobes of the impulse response and the correction factor B, for main lobes of 
the impulse response. This results in 



(0.361) 2 



,_ 7 (k-j) 2 *F (0) \ k ^ 2 



((k-i) 2 -(0.361) J )*F M (0.361) y 



20 and 



(0.361) 2 



1 — 



1-F M (0.361) 



(l-4*(0-361) 2 )*F M (0.361)' 



F M (0) 



-1 



where the normalisation factor F M (0.361) is found by substituting d = 0.361 into the 
25 equation 



kk=2 v y 
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The second approximation provides a more accurate approximation at the 
expense of a number of additional multiplication operations. 

Figure 5 gives the general schematic for a ninth-order filter with 10 filter coeffi- 
5 cients. The section surrounded by a dotted line relates to the correction factors used for 
the second approximation. This form is suitable to be implemented, for example, in a 
software program for a digital signal processor. 

A topology for both a software implementation and for the design of an inte- 
grated circuit can be found by rearranging the equations for the filter coefficients, so 
10 that a minimum number of multiplication operations using the fractional delay d are 
required. 

The equations for the filter coefficients in the first approximation (B u = 0) for 
the filter from Figure 5 are then given by a third-order polynomial in d: 



15 





*d-i(C, 


-4)*d 2 ±(C,-4)*d 3 






+ jd 2 


Td 3 


K 
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±d 3 


te 




+ id' 


+ d 3 


te 




-|d> 


±d 3 


K 
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The schematic for a design of an integrated circuit which comprises delay ele- 
ments, multiplication elements, addition elements and inversion elements, for the first 
approximation (third-order polynomial), is shown in Figure 6. The maximum approxi- 
mation error in the impulse response of this finite impulse response filter with frac- 
25 tional delay is about 4-10' 3 . 

The equations for the filter coefficients in the second approximation (B k * 0) for 
the filter from Figure 5 are then given by a fifth-order polynomial in d: 
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(h,i=iC 1 HFiC,*d-i(C 1 -4)*0+iB I )*d I ±(C 1 -4)*(l+iB,)*d 3 



+£(C, -4)*B, *d 4 +(C, -4)*B, *d 5 



5 



(h 2 ) 5 =(-l±id + |(l + iB 2 )*d z +(l + iB 2 )*d 3 -}B 2 *d 4 ±B 2 *d 5 )*C 2 
(h 3 ) 5 =( + i+id-|(l + iB 3 )*d 2 ±(l+|Bj)*d 3 + |B J *d 4 TB 3 *d 5 )*C 3 
(h 4 ) 5 =(-l±|d + |0+iB 4 )*d J +(l+{BJ*d 3 -iB 4 *d 4 ±B 4 *d 5 )*C 4 
W = (+1* {d-f (l + |B 5 )*d 2 ±(l + iB 5 )*d 3 + f B s *d< TB 5 *d 5 )*C 5 



The schematic for a design of an integrated circuit for the second approximation 
(fifth-order polynomial) is shown in Figure 7. The maximum approximation error in the 
impulse response for this finite impulse response filter with fractional delay is about 
10 2-10" 5 . 

From the formulae for the filter coefficients using the third-order polynomial 
approximation, filter coefficients can be derived for a second-order polynomial 
approximation. The filter coefficients for the side lobes of the filler using a second- 
order polynomial approximation (where subscript 2 is used to indicate the second-order 
1 5 polynomial approximation) are given by the formulae 



f h (k,d)l =(-l)"- 1 *(Kk-i)^id + (±i-(k-i))*<l 2 )*C k 



for 0 < d < 'A, and 



20 



(.MM)] =(-l) k - , *(i(k-i)T}d + (+i-(k-i))*d J )*C k 

V.i;M-i /2 



for - Vi < d < 0. For the main lobes, the filter coefficients are given by 



25 



(h k „(d,M)) I =iC 1 +(iC 1 -i)d + (2+3 + (-|±i)*C ( )*d 2 



for 0 £ d <, Vi y and 
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for - Vi < d < 0. The maximum approximation error in the impulse response is less than 
610 :j (in the equiripple embodiment), which is virtually no greater than in the case of 
the third-order polynomial approximation. 
5 Figure 8 shows a schematic for a design of an integrated circuit of the filter 

shown in Figure 5 according to a second-order polynomial approximation, derived 
from the embodiment using a third-order polynomial approximation. The equations for 
the filter coefficients in this approximation for the filter from Figure 8 are then given by 
a second-order polynomial in d: 
1 0 for positive d: 



15 



20 





-i)*d 


+ (-(W)*C 1 +2T3)*d 2 




*d 


+ (^)*d 2 )*c 2 




*d 


-(i+i)*d 2 )*C 3 




*d 


+ ( z+i)*d 2 )*C 4 




*d 


_(| T i)*d 2 >C, 


and for negative d: 






(h,l= iC,T(}C, 


-i)*d 


+ (-(i±i)*C, + 2±3)*d 2 




*d 


+(i±})*d 2 




*d 


-(f±j)*d 2 )*c 3 




*d 


+ (f±|)*d : )*C« 




*d 


-(f±i)*d 2 • )*C 5 
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The schematic shown in Figures 6 3 7 and 8 can serve as a basis both for soft- 
ware implementations of the filters and for the design of integrated circuits. 

In a further embodiment of the present invention, the inverse of the finite 
response filtering with fractional delay d is used for determining the values of the frac- 
30 tional delay d, at which the original signal has zero crossings and/or extremes. In a 
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fractional delay filter, the fractional delay parameter d (within the total phase delay) 
determines to what extent input samples are shifted in time by interpolation. For each 
clock pulse, a polynomial function in d is calculated. If the design of the fractional 
delay filter is such that the polynomial coefficients are provided directly, this filter is 
5 suitable to serve as an inverse fractional delay filter. The polynomial function and its 
derivative then directly enable the positions in time of zero crossings (d«ro) and/or of 
extremes (d ex tr) to be calculated. If the delay found at an extreme (d extr ) is introduced'as 
the input of a fractional delay filter, this directly gives the corresponding amplitude 
value as the output signal of the filter. Zero crossings can only occur between two suc- 
1 0 cessive sample values of opposite signs. To find the zero crossings it is first necessary, 
therefore, to determine where in the series of sample values this condition is met 

The exemplary schematics of Figures. 6, 7 and 8 have been elaborated in such a 
way that the separate polynomial coefficients ao, a u a 2 ... are calculated directly. From 
these coefficients, the values of the fractional delay at zero crossings can be determined 
1 5 directly and without iterations via the inverse fractional delay filter concept, and after 
differentiation, extremes of the signal such as maxima and minima can be determined. 
The calculations of these values of the fractional delay are the simpler, the lower the 
order of the polynomial approximation. The above-discussed second-order polynomial 
approximation is consequently the most suitable for use in the inverse filter with frac- 
20 tional delay. Zero crossings can therefore be found by solving the equation 

f 2 (d) = a 0 +a,*d+aj*d 2 =0 
The extremes can be found by solving the equation 

25 

f 2 '(d) = a,+2*a,*d = 0 



where the type of extreme can be found by a calculation of the second derivative if 
f 2 "(d)=2*a 2 , the extreme relating to a minimum if f 2 "(d)>0, a maximum if f 2 "(d)<0 
30 and a point of inflection if f 2 "(d)=0. 

Such determinations of zero crossings can also be carried oul with the aid of the 
third-order polynomial approximation, as discussed with reference lo Figure 6. Figure 9 
shows the schematic for a design of an integrated circuit for determining the times of 
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zero crossings and extremes and the corresponding amplitudes with the aid of the 
second-order polynomial approximation shown in Figure 8. This embodiment is 
provided with computation means 10 which perform the above-mentioned calculations 
of the zero crossings and extremes. If the fractional delay at which an extreme occurs 
(d txtr ) is directly fed back to the filter, the amplitude value at the extreme is supplied at 
the output of the filter. 

Hereinafter a number of practical applications will be discussed for using the 

method or device according to the invention. 

Hilbert transformers in the ideal case are based on infinite sine functions. This 
opens the possibility of applying the concept of filtering with fractional delay according 
to the invention, resulting in a finite impulse response by using the windowing tech- 
nique with predetermined window weights. The ideal Hilbert transformer consists of 
two filters having a flat frequency response with output signals phase-shifted by pre- 
cisely 90 degrees for common real input signals. The two output signals form a single 
complex output signal. The windowing technique according to the invention cannot be 
applied directly to the two impulse responses of the ideal Hilbert transformer which are 
currently found in the literature, as there is an odd number of samples of the sine im- 
pulse response within the even number of window weights which has to be used. This 
can be solved by causing the impulse responses of the two filters in the Hilbert trans- 
former to be mirror-antisymmetric and of the same odd order. The filter coefficients of 
the two filters are then found by substituting the value d = % (resulting in the real im- 
pulse response) and d = -V* (resulting in the imaginary impulse response) into the inter- 
polation formulae according to the invention. Real bandwidth-limited signals (for 
example coming from an antenna) are thereby directly converted into complex band- 
width-limited signals, half-band filtering being used in the process. Both time errors (or 
phase errors) and the spurious signals introduced thereby, which occur in conventional 
Hilbert transformers, are thus eliminated. 

A further application can be found for numerical controlled oscillators (NCO) 
like, for example, those used in software programs. Known NCOs, such as, for 
example, the type HSP45116 from Harris Corporation, consist of two parts, namely a 
digital K-bit phase accumulator linked to an electronic table which comprises stored 
sample values of a sine wave signal. The table contains 2 20 (about one million) numbers 
(for K=20) having a word size of 16 bits, (i.e. a total of 16 Mbit of data). By using the 
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phase accumulator it is possible to select any frequency value in the range from 0 to fJ2 
atstepsof2- 20 *f s . 

The embodiments of the fractional delay filters, as discussed in the present 
application, all have the characteristic that they exhibit very low distortion both in 
5 phase delay and in the level of the frequency response for signal frequencies close to 
zero. As a result, these are eminently suitable for constructing all sinusoidal signal 
values by interpolation of a relatively small number of sample values of a constant sine 
wave which, for example, is oversampled by a factor of 16 (b=16). Starting from a set 
of at least 2*b=32 stored sample values of a sine wave (amply sufficient for charac- 
10 terising the entire period), it is possible to construct a numerical controlled oscillator 
which is able to generate a very accurate digital sine wave of any desired frequency. 
Computer simulations have shown that the best results are achieved by using a third- 
order polynomial approximation (P0=3) and a filter order of M = 5. In this case, only 
32 sine signal samples need be stored for a better quality, compared with the known 
15 NCO (error of 1(T 7 instead of error of 1.5*10" 5 ). If use is made of filters employing a 
lower-order polynomial approximation, more sine signal samples have to stored for 
comparable quality. Simulations have shown that about 250 sine signal samples should 
be stored for a second-order (P0=2) approximation and about 2000 with a first-order 
(P0=1) approximation. 

20 If the approach is based on an interpolator with fractional delay using a third- 

order polynomial approximation, which is fed a periodic sine wave which, for example, 
is oversampled by a factor of 16 (b= 16), a set of four polynomial coefficients {a 0 , ai, 
a 2 , a 3 } is calculated at each sampling time. Owing to the periodicity of the digitised 
sine wave, not only the 32 samples of the sine wave but also the 32 related groups of 

25 four polynomial coefficients will be repeated periodically. This implies that instead of 
storing 32 sine wave samples it is better to store the 32 related groups of polynomial 
coefficients. As a result, a large number of operations of the filtering with fractional 
delay can be dispensed with. Only the calculation of the polynomial value at each 
sampling time remains. A possible set up for a numerical controlled oscillator is shown 

30 in Figure 10. The NCO comprises a K-bit phase accumulator 20 and an addressable 
table 22, for example in the form of a read-only memory (ROM). The K-bit phase 
accumulator 20 converts the desired frequency F into an addressing word and a delay 
word. The way in which this is done is assumed to be known to those skilled in the art 
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and will not be explained here in more detail. The addressable table 22 comprises 2*b 
sets of P0+1 polynomial coefficients, which are addressed by the addressing word 
which comprises the 2 log(2*b) most significant bits of the K-bit phase accumulator 20. 
The addressed P0+1 polynomial coefficients (in the examples shown a 0 , at, a 2 , a 3 ) are 
then used for a POth order approximation of a filter with a fractional delay d, formed by 
P0 multiplication elements and P0 addition elements. The fractional delay d is 
determined by the delay word which is formed by the K- 2 log(2*b) least significant bits 
of the K-bit accumulator 20. Clearly, in this embodiment the ROM 22 can be much 
smaller than in the known NCO. In the known NCO, about one million numbers need 
to be stored. If use is made of a third-order (P0=3) polynomial approximation of a fifth 
order (M=5) filter with fractional delay, only 128 polynomial coefficients need to be 
stored to achieve comparable quality. The inverter 24 in the schematic shown in 
Figure 10 changes the binary representation of the fractional delay d to a two-comple- 
ment representation. The multiplexer 26 shown in Figure 10 has been added to allow 
the choice between generating a sine or a cosine signal by feeding the two most signifi- 
cant bits of the addressing word directly to the ROM 22, or via a second inverter 28 and 
an EXOR-gate 30. The frequency of the sine/cosine is defined, in the example shown, 
by the relationship: 

frequency = (F/2 K )*f , where 0 F 2 Kl -l. 

In the case of a signal oversampled by a factor of 16, a 5th order filter with fractional 
delay and a third-order polynomial approximation (b=16, M=5, P0=3 implying four 
polynomial coefficients), the polynomial coefficients are then given by: 

a 0 (n) = -H5C 3 (x B + x._ 5 )-3C 2 (x a _, +x B _ 4 )+C ) (x n _ J +x 0 _,)} 
a,(n) = -4a 0 + 2(x n _j+x B _ 3 ) 

a,(n) = i{C 3 (-x n + x B _ 5 )+C 2 (x B _, +x B _ 4 )+C 1 (-x 0 _ 2 + x„_ 3 )} 
a 3 (n) = -4a,+4(x n . J -x 11 .3) 

where x„ (not delayed) to x„. 5 (delayed by 5 clock pulses) inclusive form a set of six 
successive sample values from a set of 2*b sine wave sample values which are calcu- 
lated according to 
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x, = sin((n + *±) * * ) n = {0,1,2,. . .,(2 * b - 1)} 

It will be obvious that instead of periodic digitised and oversampled sine wave 
shapes it is also possible to use other arbitrarily chosen periodic wave shapes which are 
limited in bandwidth and have been oversampled. 

Many more fields than the two for which examples have been elaborated here- 
inabove lend themselves to the application of the method for filtering with a fractional 
delay according to the present invention, including applications requiring time-continu- 
ous adjustable delay lines, synchronisation of receivers by using inverse filtering with 
fractional delay, etc. All these elaborations of the general concept of filtering with frac- 
tional delay are deemed to be within the scope of the appended claims. 
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Claims 



10 



1. Method for filtering a digital signal with the aid of an Mth order finite impulse 
response filtering with a fractional delay d which is a fraction of the sampling period 
associated with the digital signal, in order to obtain an interpolated value, comprising 
the step of predetermining (M+l)/2 window weights w k of a window function which is 
applied to the infinite impulse response of the ideal interpolator, characterised in that 
the window function is a symmetric stepped window function and in that the method 
further comprises the following steps: 

- determining simplified filter coefficients (hj) w o as a function of the fractional delay d 
according to the formula 



( h(M)l 



= H)"" 1 *■ 



W L 



wO 



(k-i)+d 



f h (M)^ 

j»M-i 



= (-D 



k-1 *_ 



wO 



(k-i)-d 



k = {*fW,i} = ^-i 

|d|<i 
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and subsequently applying the finite impulse response filtering to the digital signal; 
- multiplying the weighted result of the finite impulse response filtering by a common 
normalisation factor (F(d))wo. given by the formula: 



20 



(F(d)L = 



M+l 



which is a function of the simplified filter coefficients, in order lo obtain the inter- 
polated signal value from the digital signal with fractional delay d. 



25 2. Method according to Claim 1, characterised in that the window function is 
optimised for equiripple filtering, the window weights Wk being given by the formula: 



QQIV17niA1 I » 
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(2k-l) 
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3; Method according to Claim 1, characterised in that the window function is 
optimised for Lagrange filtering, the window weights w k being given by the formula: 



(\y ) = - 2 k = {1,2,...,——} 



4. Method according to Claim 1, 2 or 3, characterised in that the method is also 
suitable for |d| = l A by adding an additional common factor (Vi-d 2 ) when the simplified 
filter coefficients (hj)wo and the common normalisation factor (FCd))^ are determined. 

5. Method according to any one of the preceding claims, characterise d in that two 
finite impulse response filtering operations are employed to effect filtering which 
allows all frequencies to pass, simplified filter coefficients being determined solely 
from the window weights and the desired fractional delay d. 

6. Method according to any one of the preceding claims, characterised in that in 
determining the filter coefficients and the common normalisation factor the division is 
replaced by an approximation with the aid of one or more multiplications. 

7. Method according to Claim 6, characterised in that the approximation for the 
main lobe of the impulse response and the approximation for the side lobes of the im- 
pulse response take place separately, so that the resulting impulse response retains the 
exact values for the extreme values of the fractional delay. 

8. Method according to any one of the preceding claims, characterised in that the 
inverse of the finite impulse response filtering with fractional delay d is used for 
determining the values of the fractional delay d, at which the original signal has zero 
crossings and/or extremes. 
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9. Device for implementing the method according to any one of Claims 1 to 7 in- 
clusive, comprising delay elements, addition elements, multiplication elements and 
inverter elements for filtering a digital signal with a fractional delay d which is a frac- 
tion of the sampling period associated with the digital signal, 
5 predetermined window weights w k of a window function which is applied to the infinite 
impulse response of an ideal interpolator being implemented in multiplication elements, 
simplified filter coefficients (hj)wo being implemented as a combination of addition 
elements, multiplication elements hnd inverter elements and 

a common normalisation factor (F(d))wo being implemented as a multiplication ele- 
10 ment, 

in order to obtain the interpolated signal value from the digital signal with fractional 
delay d. 

1 0. Device according to Claim 9, characterised in that the device is further provided 
1 5 with computing means in which the inverse of the finite impulse response filtering with 
fractional delay d is employed to determine the values of the fractional delay d at 
which the original signal has zero crossings and/or extremes. 



11 . Numerically controlled oscillator for generating digitised band-limited periodic 
20 signals, comprising a K-bit phase accumulator and an addressable table characterised in 
that the addressable table comprises twice the oversampling (=2*b) of groups of P0+1 
polynomial coefficients which are addressed by the 2 log(2*b) most significant bits of 
the K-bit phase accumulator, 

the addressed P0+1 polynomial coefficients being used for a POth order approximation 
25 of a filter with a fractional delay d, formed by P0 multiplication elements and P0 addi- 
tion elements, the fractional delay d being determined by the K- 2 Iog(2*b) least signifi- 
cant bits of the K-bit phase accumulator. 



***** 
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